install.packages("strucchange")

library(strucchange)


#setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Appendix/A_Fig_6_7")





data_q<-read.csv(file="data_fig_A6_7.csv")


years<-as.numeric(as.character(unique(data_q$Year)))


num_states<-tapply(data_q$ID,data_q$Year,length)


########################### Figure 6 #################################################

pdf(file="Figure_A6.pdf")

plot(years[2:139],diff(num_states),type="b",ylab="Change in Number of States", xlab="Year")

abline(h=0,lty=2,col="red")


dev.off()

bre<-c(rep(0,80),rep(1,68))


df<-data.frame(num_states[2:139], num_states[1:138] ,bre[1:138])

names(df)<-c("num_states","lag_num_states","bre")



ft_lagdiff<-lm(diff(num_states)~lag_num_states[2:138]-1,data=df)



mean_resid_1<-round(mean(resid(ft_lagdiff)),3)

sd_resid_1<-round(sd(resid(ft_lagdiff)),3)


##Point Estimate Lammbda ##

lambda<-round(coef(ft_lagdiff),3)


##SE Point Estimate Lammbda ##

se_lambda<-round(sqrt(vcov(ft_lagdiff)[1,1]),3)



########### Print Results  ##################

out1<-data.frame(lambda,se_lambda,mean_resid_1,sd_resid_1)

names(out1)<-c("lambda","se_lambda","average residual","standard deviation of residual")

row.names(out1)<-"estimates"


print("This is the result from the model Delta Y_t = lambdaY_{t-1} + e_t")

print(out1)


###########################################################


### Note these two produce the same result ####

ft_lag<-lm(num_states~lag_num_states-1,data=df)


##Point Estimate Lammbda ##

#round(coef(ft_lagdiff),3)


##SE Point Estimate Lammbda ##


#round(sqrt(vcov(ft_lagdiff)[1,1]),3)

######################################################################################


dfb0<-df[bre==0,]

ft_lagbreak0<-lm(diff(num_states)~lag_num_states[2:nrow(dfb0)]-1,data=dfb0)

lambda_pre<-round(coef(ft_lagbreak0),3)

se_lambda_pre<-round(sqrt(vcov(ft_lagbreak0)[1,1]),3)

mean_resid_pre<-round(mean(resid(ft_lagbreak0)),3)

sd_resid_pre<-round(sd(resid(ft_lagbreak0)),3)


out2<-data.frame(lambda_pre,se_lambda_pre,mean_resid_pre,sd_resid_pre)

names(out2)<-c("lambda","se_lambda","average residual","standard deviation of residual")

row.names(out2)<-"estimates"


print("This is the result from the model Delta Y_t = lambdaY_{t-1} + e_t  - for years pre 1500")

print(out2)


######################################################################################



dfb1<-df[bre==1,]

ft_lagbreak1<-lm(diff(num_states)~lag_num_states[2:nrow(dfb1)]-1,data=dfb1)


lambda_post<-round(coef(ft_lagbreak1),3)

se_lambda_post<-round(sqrt(vcov(ft_lagbreak1)[1,1]),3)

mean_resid_post<-round(mean(resid(ft_lagbreak1)),3)

sd_resid_post<-round(sd(resid(ft_lagbreak1)),3)




out3<-data.frame(lambda_post,se_lambda_post,mean_resid_post,sd_resid_post)

names(out3)<-c("lambda","se_lambda","average residual","standard deviation of residual")

row.names(out3)<-"estimates"


print("This is the result from the model Delta Y_t = lambdaY_{t-1} + e_t  - for years post 1500")

print(out3)




########## Figure 7 #########



pdf(file="Figure_A7.pdf")



f_out<-Fstats(diff(num_states)~lag_num_states[2:138]-1,data=df,from=.05,to=.9)

plot(f_out)


dev.off()

